A STABILIZED FINITE ELEMENT METHOD FOR 
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Abstract. A recently developed Eulerian finite element method is applied to solve advection- 
difl^usion equations posed on hypersurfaces. When transport processes on a surface dominate over 
difl'usion, finite element methods tend to be unstable unless the mesh is sufficiently fine. The paper 
introduces a stabilized finite element formulation based on the SUPG technique. An error analysis of 
the method is given. Results of numerical experiments are presented that illustrate the performance 
of the stabilized method. 
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1. Introduction. Mathematical models involving partial differential equations 
posed on hypersurfaces occur in many applications. Often surface equations are cou- 
pled with other equations that are formulated in a (fixed) domain which contains the 
surface. This happens, for example, in common models of multiphase fluids dynamics 
if one takes so-called surface active agents into account [T]. The surface transport 
of such surfactants is typically driven by convection and surface diffusion and the 
relative strength of these two is measured by the dimcnsionlcss surface Peclet number 
Pcs = jj^- Here U and L denote typical velocity and lenght scales, respectively, and 
Ds is the surface diffusion coefficient. Typical surfactants have surface diffusion coef- 
ficients in the range Dg ~ 10"'^ — 10~^ cm? / s [2], leading to (very) large surface Peclet 
numbers in many applications. Hence, such applications result in advection-diffusion 
equations on the surface with dominating advection terms. The surface may evolve 
in time and be available only implicitly (for example, as a zero level of a level set 
function) . 

It is well known that finite element discretization methods for advection-diffusion 
problems need an additional stabilization mechanism, unless the mesh size is suffi- 
ciently small to resolve boundary and internal layers in the solution of the differential 
equation. For the planar case, this topic has been extensively studied in the literature 
and a variety of stabilization methods has been developed, see, e.g., [3]. We arc, how- 
ever, not aware of any studies of stable finite clement methods for advection-diffusion 
equations posed on surfaces. 

In the past decade the study of numerical methods for PDEs on surfaces has been 
a rapidly growing research area. The development of finite element methods for solv- 
ing elliptic equations on surfaces can be traced back to the paper [3] , which considers 
a piecewise polygonal surface and uses a finite element space on a triangulation of this 
discrete surface. This approach has been further analyzed and extended in several 
directions, see, e.g., [S| and the references therein. Another approach has been intro- 
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duced in [B] and builds on the ideas of [3 • The method in that paper apphcs to cases 
in which the surface is given imphcitly by some level set function and the key idea is to 
solve the partial differential equation on a narrow band around the surface. Unfitted 
finite element spaces on this narrow band arc used for discretization. Another surface 
finite element method based on an outer (bulk) mesh has been introduced in [8] and 
further studied in [9l [10] . The main idea of this method is to use finite element spaces 
that are induced by triangulations of an outer domain to discretizc the partial differ- 
ential equation on the surface by considering traces of the bulk finite element space 
on the surface, instead of extending the PDE off the surface, as in [7J [B] • The method 
is particularly suitable for problems in which the surface is given implicitly by a level 
set or VOF function and in which there is a coupling with a differential equation in 
a fixed outer domain. If in such problems one uses finite clement techniques for the 
discetization of equations in the outer domain, this setting immediately results in an 
easy to implement discretization method for the surface equation. The approach does 
not require additional surface elements. 

In this paper we reconsider the volume mesh finite element method from [S] and 
study a new aspect, that has not been studied in the literature so far, namely the 
stabilization of advection-dominated problems. We restrict ourselves to the case of a 
stationary surface. To stabilize the discrete problem for the case of large mesh Peclet 
numbers, we introduce a surface variant of the SUPG method. For a class of stationary 
advcction-diffusion equation an error analysis is presented. Although the convergence 
of the method is studied using a SUPG norm similar to the planar case [3], the analysis 
is not standard and contains new ingredients: Some new approximation properties for 
the traces of finite elements are needed and geometric errors require special control. 
The main theoretical result is given in Theorem 13.101 It yields an error estimate in 
the SUPG norm which is almost robust in the sense that the dependence on the Peclet 
number is mild. This dependence is due to some insufficiently controlled geometric 
errors, as will be explained in section [3.71 

The remainder of the paper is organized as follows. In section [51 we recall equa- 
tions for transport-diffusion processes on surfaces and present the stabilized finite 
element method. Section [3] contains the theoretical results of the paper concern- 
ing the approximation properties of the finite element space and discretization error 
bounds for the finite element method. Finally, in section[4|results of numerical experi- 
ments are given for both stationary and time-dependent advection-dominated surface 
transport-diffusion equations, which show that the stabilization performs well and 
that numerical results are consistent with what is expected from the SUPG method 
in the planar case. 

2. Advection-difFusion equations on surfaces. Let be an open domain in 
M'^ and F be a connected compact hyper-surface contained in O. For a sufficiently 
smooth function g ; f2 — > R the tangential derivative at F is defined by 



where nr denotes the unit normal to F. Denote by Ar the Laplace-Beltrami operator 
on F. Let w : — > R"^ be a given divergence-free (div w = 0) velocity field in O. If the 
surface F evolves with a normal velocity of w • nr, then the conservation of a scalar 
quantity u with a diffusive fiux on F(t) leads to the surface PDE: 



Vr.g = - (V.g • nr)nr, 



(2.1) 



u + (divr w)u — sAyu = 



on F(i), 



(2.2) 
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where u denotes the advective material derivative, e is the diffusion coefficient. In 
[TTj the problem (|2.2I) was shown to be well-posed in a suitable weak sense. 

In this paper, we study a finite element method for an advection-dominated prob- 
lem on a steady surface. Therefore, we assume w • nr = 0, i.e. the advection velocity 
is everywhere tangential to the surface. This and div w = implies divr w = 0, and 
the surface advection-diffusion equation takes the form: 

U4 -I- w • Vrw - eArw = on T. (2.3) 

Although the methodology and numerical examples of the paper are applied to the 
equations (|2.3p . the error analysis will be presented for the stationary problem 

- eArw + w • Vpu + c(x)u = / on T, (2.4) 

with / e i^(r) and c(x) > 0. To simplify the presentation we assume c(-) to be 
constant, i.e. c{x) — c>Q. The analysis, however, also applies to non-constant c, cf. 
section [3771 Note that ()2.3|) and (j2.4|) can be written in intrinsic surface quantities, 
since w • Vpu = vvr • Vpit, with the tangential velocity wr = w — (w • nr)nr. 
We assume wr € H^'°°{T) n L°°iT) and scale the equation so that || v^fr||L°=(r) = 1 
holds. Furthermore, since we are interested in the advection-dominated case we take 
£ E (0, 1]. Introduce the bilinear form and the functional: 

a{u,v) Vrit-Vrwds + J {w-\7ru)vds + J cuwds, 

f{v) := J fvds. 



The weak formulation of (|2.4I) is as follows: Find u £ V such that 

a{u,v) = f{v) yveV, (2.5) 



with 

V 



{u e iJi(r) I /pwds = 0} ifc = 0, 

iji(r) ifoo. 



Due to the Lax-Milgram lemma, there exists a unique solution of (j2.5p . For the case 
c = the following Friedrich's inequality [12] holds: 

||«lli2(r) < ||Vr«||i.(r) V € K (2.6) 

2.1. The stabilized volume mesh FEM. In this section, we recall the volume 
mesh FEM introduced in [8] and describe its SUPG type stabilization. 

Let {Th}h>o be a family of tetrahedral triangulations of the domain f2. These 
triangulations are assumed to be regular, consistent and stable. To simplify the 
presentation we assume that this family of triangulations is quasi-uniform. The latter 
assumption, however, is not essential for our analysis. We assume that for each Th a 
polygonal approximation of F, denoted by F/i, is given: Th is a C*^'^ surface without 
boundary and Th can be partitioned in planar triangular segments. It is important to 
note that F^ is not a "triangulation of F" in the usual sense (an 0{h^) approximation 
of F, consisting of regular triangles). Instead, we (only) assume that Th is consistent 
with the outer triangulation Th in the following sense. For any tetrahedron St S Th 
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such that nicas2(5'T nr^) > 0, define T — STC\Th. We assume that every T e F/j is a 
planar segment and thus it is either a triangle or a quadrilateral. Each quadrilateral 
segment can be divided into two triangles, so we may assume that every T is a triangle. 
An illustration of such a triangulation is given in Figure 12.11 The results shown in 
this figure are obtained by representing a sphere F implicitly by its signed distance 
function, constructing the piecewise linear nodal interpolation of this distance function 
on a uniform tetrahedral triangulation 7/i of f2 and then considering the zero level of 
this interpolant. 



Fig. 2.1. Approximate interface for a sphere, resulting from a coarse tetrahedral triangula- 
tion (left) and after one refinement (right). 

Let Th be the set of all triangular segments T, then Th can be decomposed as 

F,, = y T. (2.7) 

Note that the triangulation J-h is not necessarily regular, i.e. elements from T may 
have very small internal angles and the size of neighboring triangles can vary strongly, 
cf. Figure [Ol In applications with level set functions (that represent F implicitly), 
the approximation F;, can be obtained as the zero level of a piecewise linear finite 
element approximation of the level set function on the tetrahedral triangulation Th ■ 

The surface finite element space is the space of traces on Th of all piecewise linear 
continuous functions with respect to the outer triangulation Th ■ This can be formally 
defined as follows. We define a subdomain that contains F^: 

oJh= U St, (2.8) 

an a corresponding volume mesh finite element space 

Vh := {vh e CiuJh) I Vh\s^ € Pi V T e Th}, (2.9) 

where Pi is the space of polynomials of degree one. Vh induces the following space on 
Th-. 

V[ := {^Jh e H^Vh) \3vh£ Vh such that = «h|rj. (2.10) 

When c = 0, we require that any function Vh from V^ satisfies Jp^ Vh ds = 0. Given 
the surface finite clement space V^ , the finite element discretization of (|2.5p is as 
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follows: Find Uh G such that 

£ / Vr;,M • Vr;,wds + / (w" • Vr;,M)wds + / cuvds = / fvds (2.11) 

JTh JTh JTh JPh 

for all Vh G . Here w*^ and arc the extensions of wr and /, respectively, along 
normals to F (the precise definition is given in the next section). Similar to the 
plain Galerkin finite element for advection-diffusion equations, the method (|2.1ip is 
unstable unless the mesh is sufficiently fine such that the mesh Peclet number is less 
than one. 

We introduce the following stabilized finite element method based on the standard 
SUPG approach, cf. [5]: Find Uh € such that 

ahiuh.Vh) = fh{vh) V Vh e V^, (2.12) 

with 

aji{u,v) :=e / Vr,,w • Vr,,u ds + / cuvds 

+ \ I (w'= • Vr,.u)wds- / (w*^ • Vr^t;)uds (2.13) 
+ ^ 5t I (— eAr^u + w*^ • Vr^u + cu)w'^ • Vr^wds, 

Mv):^ f rvds+ I r(w^.Vr,i;)ds. (2.14) 

The stabilization parameter St depends on T C St- The diameter of the tetrahedron 

'iStII"*^'^II-L°°(T) 

St is denoted by /i^^. Let Pe^ := ^ — be the cell Peclet number. We 

take 

il Pej- > 1, 



and (5t = min{5T,c-^}, (2.15) 



if Per < 1, 



with some given positive constants (5o, > 0. 

Since £ is linear on every T we have Ap^w/i = on T, and thus ah{uh,Vh) 
simplifies to 

ahiuh,Vh) = s Vr,,ith-Vr^w/ids + i / (v^f^ • Vr,,it/i)w/,. - (w" • Vr,,u;i)uft ds 
+ f cuh{vh+Si^)w^ ■Vr^vh)ds+ f 5(x)(w^ • Vr,M^)(W- • V^,^;/^) ds, (2.16) 

where 5{x.) = 6t for x S T. 

3. Error analysis. The analysis in this section is organized as follows. First we 
collect some definitions and useful results in section [3Tl In section [321 derive a 
coercivity result. In section 15751 we present interpolation error bounds. In sections l3.4l 
and l3.5[ continuity and consistency results are derived. Combining these analysis we 
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obtain the finite element error bound given in section 13.61 In the error analysis we 
use the foUowing mesh-dcpendcnt norm: 

\\u\\,:=(ef |Vr,u|2ds+ / J(x)|w^ • Vr.wp ds + / cjupdsV. (3.1) 

Here and in the remainder | • | denotes the Euclidean norm for vectors and the corre- 
sponding spectral norm for matrices. 

3.1. Preliminaries. For the hypcrsurface F, we define its /i- neighborhood: 

Uh ■■= {x G I dist(x, F) < coh}, (3.2) 

and assume that cq is sufBciently large such that ojh C Uh and h sufficiently small 
such that 

Scoft. < (maxj=i,2||Ki||Loo(r)) ^ (3.3) 

holds, with Ki being the principal curvatures of F. Here and in what follows h denotes 
the maximum diameter for tctrahcdra of outer triangulation: h = max diam(S'). 

Let d : [//i — ?► M be the signed distance function, \d{x.)\ = dist(x, F) for all x £ Uh- 
Thus F is the zero level set of d. We assume c? < in the interior of F and d > Q 
in the exterior and define n(x) \/d{x) for all x e Uh- Hence, n = np on F and 
|n(x)| = 1 for all x e Uh- The Hessian of d is denoted by 

H(x) := W^d{x) e R3''^ X e Uh- (3.4) 

The eigenvalues of H(x) are denoted by ki(x), K2(x), and 0. For x e F the eigenvalues 
Ki,i = 1,2, are the principal curvatures. 

For each x e Uh, define the projection p : [//i — > F by 

p(x) =x-rf(x)n(x). (3.5) 

Due to the assumption p.3p . the decomposition x = p(x) + (i(x)n(x) is unique. We 
will need the orthogonal projector 

P(x) ;= I - n(x)n(x)^, for x € Uh- 

Note that n(x) = n(p(x)) and P(x) = P(p(x)) for x G Uh holds. The tangential 
derivative can be written as Vrg(x) ~ PV5(x) for x e F. One can verify that for 
this projection and for the Hessian H the relation HP = PH = H holds. Similarly, 
define 

P/i(x) := I — nr,, (x)nr^(x)^, for x G Th, x is not on an edge, (3.6) 

where nr^ is the unit (outward pointing) normal at x G Th (not on an edge). The 
tangential derivative along Th is given by Vr,,g(x) = P/i(x)V.g(x) (not on an edge). 

Assumption 3.1. In this paper, we assume that for ah T e Th- 

ess supxg-r|d(x)| < ci/i|^, (3.7) 
ess supxgr|n(x) - nrjx)| < C2/ist> (3-8) 
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where hsj, denotes the diameter of the tetrahedron St that contains T, i.e., T = 
St n Fft and constants ci, C2 arc independent of h, T. 

The assumptions (|3.7p and (|3.8p describe how accurate the piecewise planar approxi- 
mation Th of r is. If F/i is constructed as the zero level of a piecewise linear interpola- 
tion of a level set function that characterizes F (as in Fig. l2.ip then these assumptions 
are fulfilled, cf. Sect. 7.3 in [J. 

In the remainder, A< B means A < cB for some positive constant c independent 
of h and of the problem parameters e and c. A B means that both A < B and 
B<A. 

For X e Tk , define 

/i^(x) = (1 - d(x)Ki(x))(l - (i(x)K2(x))n^(x)n/i(x). 

The surface measures ds and ds/i on F and Th, respectively, are related by 

/^,,(x)dsft(x) = ds(p(x)), X e Th. (3.9) 

The assumptions (13.71) and p. 81) imply that 

ess sup^grjl - ^^h)<h^, (3.10) 

cf. (3.37) in [8]. The solution of (|2.4p is defined on F, while its finite element approx- 
imation Uh G is defined on Th- We need a suitable extension of a function from F 
to its neighborhood. For a function w on F we define 

w'=(x) := w(p(x)) for all X e [/,,. (3.11) 

The following formula for this lifting function are known (cf. section 2.3 in |13)): 

Vu^(x) = (I - d(x)H)Vru(p(x)) a.e. on Uh, (3.12) 
Vr,w^(x) - P4x)(I - rf(x)H)Vru(p(x)) a.e. on F;,, (3.13) 

with H = H(x). By direct computation one derives the relation 

W^u^ix) = (P - d(x)H)V^u(p(x))(P - rf(x)H) - (n^Vru(p(x))H 

- (HVru(p(x)))n'^ - n(HVru(p(x)))^ - dVpH : Vru(p(x)). (3.14) 

For sufficiently smooth u and < 2, using this relation one obtains the estimate 

l^'^w'WI < I l^r"(PW)l + |Vrw(p(x))| I a.e. on Uh, (3.15) 
\ImI=2 / 

(cf. Lemma 3 in [3]). This further leads to (cf. Lemma 3.2 in [5]): 

WD^W'WL^iu,) < v^Mh-^^d, Ia^I < 2. (3.16) 

The next lemma is needed for the analysis in the following section. 
Lemma 3.1. The following holds: 

l|divrftW||i==(rh) ^ h\\^r'w\\L^(^r)- 
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Proof. We use the following representation for the tangential divergence: 

divr w(x) = tr(Vrw(x)) = tr(PVw(x)). (3.17) 
Take x e F/i, not lying on an edge. Using p.l2p we obtain 
divr„w''(x) 

= tr(P;,Vw'=(x)) = tr (P,,(I - rf(x)H)Vrw(p(x))) 

= tr (PVrw(p(x))) + tr ((P^. - P)Vrw(p(x))) - d(x) tr (P,,HVrw(p(x))) . 

The first term vanishes due to tr (PVrw(p(x))) = divrw(p(x)) = 0. The second 
and the third term can be bounded using p.7|) . (|3.8p : 

\Ph-P\<h, \d{^)Phll\<h\ 

This proves the lemma. □ 

3.2. Coercivity analysis. In the next lemma we present a coercivity result. 
We use the norm introduced in (|3.1I) . 

Lemma 3.2. The following holds: 

ah{vh,Vh) > ^\\vh\\l for all Vh £ . (3.18) 

Proof. For any Vh € , we have 

ah{vh,Vh) ^\\vh\\l+ C(5(x)t;,i(w'= • Vr^w/i)ds. (3.19) 

The choice of St, cf. (|2.15p . implies ci5(x) < 1. Hence the last term in p.lQp can 
be estimated as follows: 

I / cS{x.)vh{w'' ■ Vr^Wh)ds| 
This yields (jXTS)) . □ 

As a consequence of this result, we obtain the well-posedness of the discrete problem 

(12321). 

3.3. Interpolation error bounds. Let Ih ■ C{uih) Vh be the nodal interpo- 
lation operator. For any u £ i^^(F) the surface finite element function (/ftu'^)|r^ G 

is an interpolant of u'^ in . 

For any u £ i^^(F) the following estimates hold [8]: 

\\u^-ihu')\rJmr,)<h^u\\HHr), (3.20) 
W'^r.u^ -^r,Ahu')\rJ\mr,) < hWuWHHD- (3.21) 

Using these results we easily obtain an interpolation error estimate in the || • ||*- 
norm: 
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Lemma 3.3. For any u G -ff^(r) the following holds: 

U - {IiX)\tAU < He^'^ + h^'^ + cH)||u|lH2(r). (3.22) 

Proof. Define Lp u" ~ {Ihu'')\r^ G H^{Th). Using tlie definition ((^1^ of 5(x), 
we get 

/ <5(x)|w^.Vr,^|'ds= [ Srl^^-VrM^ds 
■'^^ TerjT (3.23) 

Tfie remaining two terms in Hu*^ — (/hU*^)!]* are estimated in a straiglitforward way 
using ([5:^ and Tfiis and imply tfie inequality ((X^ . □ 

The next lemma estimates the interpolation error on the edges of the surface 
triangulation. In the remainder, £h denotes the set of all edges in the interface trian- 
gulation J-h- 

Lemma 3.4. For all u e H'^{T) the following holds: 

( E X^"' - ^^"')lr. ds] < /.'/'hiU.(r). (3.24) 

Proof Define u'^ - Ihu" E H^{uJh)- Take E e £h and let T e he 
a corresponding planar segment of which E is an edge. Let be a side of the 
tetrahedron St such that E C W. From Lemma 3 in [Uj we have 

\mhiE)<h-'Mh^w) + m\miwy 
From the standard trace inequality 

Mhi0ST)<h-^MhiST)+'^MHHST) for all w € H\St), 
applied to cj) and dxi4', * = 1, 2, 3, we obtain 

hmhiw)<Mmis^)+h'\K\\hiSrr 
From standard error bounds for the nodal interpolation operator we get 

Summing over E e £h and using ||M'^||_ff2(^j^) < /i^||m||^2(p), cf. p.l6p . results in 



,3|L.||2 



which completes the proof. □ 
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Fig. 3.1. 

3.4. Continuity estimates. In this section we derive a continuity estimate for 
the bihnear form a/i(-, •). If one applies partial integration to the integrals that occur 
in ah{-,-) then jumps across the edges E € £h occur. We start with a lemma that 
yields bounds for such jump terms. Related to these jump terms we introduce the 
following notation. For each T G Th, denote by mhls the outer normal to an edge E 
in the plane which contains element T. Let [mft,]!^; = + be the jump of the 
outer normals to the edge in two neighboring elements, c.f. Figure [33] 

Lemma 3.5. The following holds: 

|P(x)[m,i](x)| < h"^ a.e. x£E. (3.25) 

Proof. Let E be the common side of two elements Ti and T2 in Th, and n^, 
, and are the unit normals as illustrated in Figure 13.11 Denote by s^, 
the unit (constant) vector along the common side E, which can be represented as 
s;i = n^ X = x n^ . The jump across E is given by 

[m,,] = Sh X (n+ - n"). 

For each x. E E and p(x) £ F, let n = n(p(x)) be the unit normal to F at p(x) and 
P = P(x) = I — nn-^ the corresponding orthogonal projection. Using p.8|) . wc get 

- n+l < |n+ - n| + |n^ - n| < hs^.^ + hs^^_ < h. 

Since |n^| = |n^| = |n| = 1, the above estimate implies 

n^ — n^ = c/i^n + ei, ei _L n, |ei| < h. 

We also have 

Sh = n+ X m+ = (n + (n+ - n)) x m+ n x + ea, jea] < h. 
We use the decomposition 

P[m,i] = P [(n X m+ + eg) x [ch^n + ei)] . 

Since ei _L n we have (n x m^) x ei || n and thus P((n x m^) x ei) = 0. Therefore, 
we get 

\P[mh]\<h^ + \ei\\e2\<h\ (3.26) 
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i.e., the result ((X^ holds. □ 

In the analysis below, we need an inequality of the form ||wh|lL2(r^) ^ for ^-H 

Vfi e . This result can be obtained as follows. First we consider the case c = 0. 
Then the functions Vh G satisfy Jp ii/i ds = 0. We assume that in a discrete 
analogon of the Friedrich's inequality ()2.6p holds uniformly with respect to h, i.e., 
there exists a constant C p independent of h such that 

hhWhir,} < C^F||Vr,i^,.ili2(r,) for all Vh G . (3.27) 

Now we reduce the parameter domain £ e (0, 1], c > as follows. For a given generic 
constant cq with < cq < 1, in the remainder we restrict to the parameter set 



e e (0, 1], c e {0} U [coe,oo). 



(3.28) 



For c > we then have 



2co ,i9 2 „, 2 

< ^-Ilc3w,||?,,,_, < -, \\v,,\\l < — 



\vh\\l 



and combing this with the result in p.27p for the case c = we get 

||i;;,|U.(r,) < -l=\\v4. for aU G (3.29) 
Ve + c 

and arbitrary e g (0, 1], c € {0} U [cqE, oo). 



In the proof of Theorem 13.81 below we need a bound for ||wh||L2(-£^) in terms of 
I It; /ill*. Such a result is derived with the help of the following lemma. 

Lemma 3.6. Assume the outer tetrahedra mesh size satisfies h < ho, with some 
sufficiently small ho ~ 1, depending only on the constant C2 from p.8|) . The following 
holds: 

J2 f vlds<h-'\\v4l,^r,)+h\\^r,v^\\hir,) forallv^eV^ (3.30) 



Proof. Let E ^ £hhe an edge of a triangle T G Th and St ^ Th is the correspond- 
ing tetrahedron of the outer triangulation. Consider the patch w(S't) of all S E Th 
touching St- Denote uj{St) ~ 2;(5T)nF;i. Let Vh be an arbitrary fixed function from 
. We shall prove the bound 

/ «^ds</^-i||i>,||i.(^(5^))+/i||Vr,z;,.||i.(^(s,)). (3.31) 
Je 

Then summing over all E € £fi and using that uj{St) consists of a uniformly bounded 
number of tetrahedra (due to the regularity of the outer mesh) , we obtain (|3.30|) . 

Let P be a plane containing T. We can define for sufficiently small h an injcctive 
mapping cj) : ijj{St) — >■ P such that |V0| < 1 and |V(^~^)| < 1. For example, </> can be 
build by the orthogonal projection on P. Then |V(/)| < 1 and |V((/)^^)(x)| < (sina)"-'^, 
where a is the angle between P and n/i(0^^(x)). Due to assumption (|3.8p we have 
1 < sin a for sufficiently small h. If (f) is the orthogonal projection on P, then (t){E) = 
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E. Thus we get 



vlds= / {vhocf) i)2ds, 

E J<p{E) 
\\vh\\L^uj{ST)) - hh o ^~^|iL2(0(c^(ST)))> 
|Vr„l'/i||L2(a;(5T)) - II W K O ) || ^2 (^(^(5.^))) . 



(3.32) 



Due to the shape regularity of 5* € uj{St) we have 

h < dist(£;,a5(5T)) < dist{E, duj (St)). 

Hence, from |V(/>^i| < 1 it follows that h < dist{4:{E),d (I){uj{St)))- Thus, we may 
consider a rectangle Q C 0(a;(5'T)) such that E = (j){E) is a side of Q and |Q| ~ h\E\. 
By the standard trace theorem and scaling argument we get 



(vho^ 'fds<h ^\\vhO(f> '\\l2iQ)+h\\Vp{vhO(l) i)||i2(Q). 



This together with ([5:^ and Q C (P{uj{St)) implies ([53T|) . □ 

An immediate consequence of the lemma and p.29p is the following corollary. 
Corollary 3.7. The following estimate holds: 



hY. I «'ds< + forallvH€V[. 



(3.33) 



E£Sh 



We are now in position to prove a continuity result for the surface finite element 
bilinear form. 

Lemma 3.8. For any u e H^{T) and G , we have 



Proof. Define (j) = u'^ — {Ihu'^)\Thi then 



(3.34) 



1 



+ / 9 ((^'^ ■ '^r,»w/i - (w^ • VT^Vh)(l>) +c(j)Vh ds 



(3.35) 



We estimate ah{(t>,Vh) term by term. Due to p.20p . p.2ip . we get for the first term 
on the righthand side of p.35p : 



£/ Vr,> • Vr^uftds < eh\\%i\\H2{T)\\"^THVh\\L-2{T„) < £''-h\W\\H'2{r)\\vh\\*- (3.36) 
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To the second term on the righthand side of p.35p we apply integration by parts: 

i((w'' • Vr»Uh - (w*" • Wr^Vh)4>) + o^w^ds 
h ^ 

-'r;, Jv^ (3.37) 

+ 7: / (w"" • m,,)0i;,ids - - / (divr,, w'=)(/)iJ,ids 

= :/l+/2+/3+/4. 

The term 1\ can be estimated by 

l-^il ^ c^|l(?!)||i2(r^)||Vcw;,.||L^(r,o ^ /»^c^ ||u||ff2(r)||v/i||*. 

To estimate I2 , we consider the adveetion-dominatcd case and the diffusion-dominated 
case separately. If Per > 1, we have 

1/2 



<max(/i-i/2,ci/2)||0j|^,(^^ (^^5j,(w^.Vr,f„)2ds^ ^ , 



and if Per < 1: 



(w'' • Vr^u/i)(/)ds < ||w*=||ioo(T)||Vr„-y/i||L2(T)ll0llL2(T) 
Summing over T G J^i we obtain 

The term /a is estimated using Pw*^ = w"^ , Lemmas 13.41 13.51 f^nd Corollary 13.71 
l^3|<| E / (W- [m^])(/.t^,.ds 

1/2 / ^ 1/2 

(3.38) 



^ ds 

1/2 



The term /4 in p.37p can be bounded due to Lemma 13. 1[ the interpolation bounds 
and dSSni): 

I/4I < ^l|divr^w'=||L»(r,,)||0|ii2(r^)||uft||i2(r^) < h^\\u\\H--i^Y)\\vh\\L^{Tn) 
< -j==M\H-^r)\\vh\\*- 
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Finally wc treat the third term on the righthand side of (|3.35l) . Using (5T||W^||ioo(-y) < 
^! S 1> ^tc < 1 and the interpolation estimates p.20[) and (|3.2ip we obtain: 



^ (5t / (-eAr;,0 + • Vr,> + c<j>)w'' ■ Vr^Vh ds 



J2 ST{e^Ar,X\\UT) + \\^''-^r,MlHT)+c^^l^T))] IM* ^^'^^^ 

< ie'/^ + h^^^ + c-2h)h\\u\\H2^r)\\vh\\*. 
Combing the inequalities p.36p - p.39p proves the result of the lemma. □ 

3.5. Consistency estimate. The consistency error of the finite element method 
(|2.12p is due to geometric errors resulting from the approximation of T by Th- To 
estimate this geometric errors we need a few additional definitions and results, which 
can be found in, for example, [T3|. For x e F^i define Ph{^) = I— n;i(x)n(x)^/ (n;i(x)- 
n(x)). One can represent the surface gradient of m G H^{T) in terms of VrhU*^ as 
follows 

Vru(p(x)) = (I - d(x)H(x))-ip,,(x)Vr,u"(x) a.e. x e F,,. (3.40) 
Due to p.gp . we get 

[vruWrvds^ A,,Vr^u"'Vr^u'' ds for all w € ^^(F), (3.41) 

with A,,(x) = /Xft(x)P^(x)(I-d(x)H(x))-2p,^(x). From w-n = onF and w^(x) = 
w(p(x)), n(x) = n(p(x)) it follows that n(x) • w'^(x) = and thus w(p(x)) = 
P/j(x)w'^(x) holds. Using this, we get the relation 

/"(w Vrit)wds = /" (BhW^ • Vr^u'=)w'=ds, (3.42) 

with Bh, = ^;,,(x)P^(/ — dH.)^^Ph. In the proof we use the lifting procedure F^, — > F 
given by 

«,Up(x)) :=v/,(x) for xeF,,. (3.43) 

It is easy to see that € _ff^(F). 

The following lemma estimates the consistency error of the finite element method 
(I2T21) . 

Lemma 3.9. Let u G H^{T) be the solution of (|2.5p , then we have 

\fh{vh) - ah{u'',Vh)\ ^ f 1 1 h - 

sup r, — r. < (fe2 +c-^h+ ^— ) fe( M g2(r) + / L^(r) ■ (3.44) 



Proof. The residual is decomposed as 

fh{vh) - ah{u^,Vh) = fh{vh) - /(w'j + a(M,W/J - ah(ii'',w/J. 



(3.45) 
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The following holds: 

fivh) = / /Wft ds = / fihfvh ds, 
a{u,vl^)—e J VrwVrwJj ds + J (w • Vrw)i;^ds + J cuv'/^ds 

= ^ j VrwVrwjids+i j (w • Vru)i4 - (w • Vrujjwds + j cuv\^ds 

- \ I (B/.w'' • Vr„i'/i)u''ds + / fihcu^vhds. 
Substituting these relations into p.45p and using (|2.13p . (|2.14p results in 

fh{vh) ~ ah{u'',Vh) ^ / {1 - ^ih).rvhds + e (A/^ - P/jVr,,-"" • Vr,j'/i ds 
+ 1 f ((B^-P^)w^.Vr,u'=Kds-i / ((B„-P,Ow^.Vr,«/.Kds 

=: h+h+h+h + h+ni. (3.46) 



We estimate the /j terms separately. Applying (|3.10p and p.29p we get 

h < /^'ll/1lL^(r.)l!"/.llL^(r.) < ^^||/||L^(r)ll«/.||*, (3.47) 

h < h^c^\\u''\\L2^rJ\Vcvh\\L^(r,,) < h'^c^\u\\L2^r)\\vh\\*. (3.48) 
One can show, cf. (3.43) in [8], the bound 

|P^- A;,| = \Ph{I-Ah)\<h\ 

Using this we obtain 

h < eh^Vr,u^U-HT,)\\^r,M\LHr,) < s'^'h'\\u^\\H^r)\\vh\U. (3.49) 
Since (I — dH)^^ = I + 0{h^), we also estimate 

|B,,-P,,| </i2 + |A,,-P,,| </l2. 

This yields 

/3 < /^'llVr.i^1|L^(r.)||i'/.||L^(r.) < -^\\u\\H^r)\\vh\U. (3.50) 

To estimate I4 we use the definition (|2.15p of 6t- If Per < 1, then e~5 |lw^|li=o(7-') < 
V 2j|w"^||2^^j,^/i^^ holds. If Pey > 1 , then ^ < max(c2 , ^ ||w'^||£oo^j,^ft,g^ ) holds. 
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Using the assumption that the outer triangulation is quasi- uniform we get h^^ ^ h^^- 
Thus, wc obtain 

min{e~^ ||w'^||i=o('r), ^ } < max(c^,/i~5) < + , 

and 

/4 <max |B/i - PftI (en|Vr;,i'/i||L2(T) +<5|.||w^^ • Vr^w/i!lL2(T)) 

X min{e ^ \\^^''■\\^^^^J.^^,SJ,^}\\u''■\\L2(^T■^ 
< h{c^ h + h^)\\vh\\4u\\L^{r)- (3.51) 
Now we estimate ili. Using the equation —eAru + w • Vrw + cw = / on F we get 

TTi = ^ (5t / (/ o p + eAr^u" - • Vr^M" - cw'")w'' • Vr^w/i ds 

= ^ j (-£(Arw)op + eAr^u'')w''- Vr;,w/,.ds 

+ ^ (5t / ((w*^ • Vru) o p - w"^ • VrhW'^)w'^ • Vr^w/i ds 

+ ^ (5t / (cit op — cm'°)w'^ • Vr^w?! ds 

=: nl + n'l + nl (3.52) 

From = It o p it foUows that TTf = holds. For Ilf we obtain, using p.40p and 
l/i^iR,, - P,,| < 1^,, - l||^,,|-i|B;,| + \Bh - Ph\ < h^, 

</tM ^ STllw^WUlVr^u^LHT)] ( E '5Tl|w^ -Vr,^'.||i.(T)) 

<hi\\u\\H2ir)\\vh\\*. (3.53) 
Since 5t£ ^ h'^, we get 



^2 < E e''5T|l(Ar^i)op-Ar,^11i.(^) E '^tIIw^ • Vr,^;,.|li.(^) 

<teH 5] ||(Ari.)op-Ar,^.'=||i.(^)) \\vy\\... (3.54) 

We finally consider the term between brackets in p.54p . Using the identity divp f = 
tr(Vrf) and n • Vu"^(p(x)) = we obtain for x G T, with := VV^. 

Aru(p(x)) = divr Vru(p(x)) = tr(PVPVu"(p(x))) = tr(PV2u'=(p(x))P). (3.55) 
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From the same arguments it follows that 

Ar,u"(x) = tr(P,,VPftVu'--(x)) = tr(Ph'^^u''{^)Ph) (3.56) 

holds. Using (|3Ta and |d(x)| < h^, |P - P,,| < h, |H| < 1, |VH| < 1 we obtain 

PhV^u%^)Fh = PV2u-(p(x))P + R, 

|R|</.(|VV(p(x))| + |V^.'=(p(x))|). 

Thus, using p.55p and (|3.56p . we get 

|Ar«(p(x)) - Ar,u^(x)| < |tr(PV2w^(p(x))P - P,V2««(x)P,)| 

</i(|W(p(x))| + |Vi.'=(p(x))|), 

and combining this with (|3.54p yields 

nl<he^ J2 ll(^rw)op- Ar,M^||i2(7,) J \\vh\\, < h^\\u\\H2(^r)\\vh\\*. 

Combining this with the resuhs (|06)) - (|33T|) and ((3?53)) proves the lemma. □ 

3.6. Main theorem. Now we put together the results derived in the previous 
sections to prove the main result of the paper. 

Theorem 3.10. Let Assumvtion \3. 1\ be satisfied. We consider problem parame- 
ters s and c as in (j3.28p . Assume that the solution u of (|2.5p has regularity u £ H^{T). 
Let Uh be the discrete solution of the SUPG finite element method (|2.12[) . Then the 
following holds: 

\\u' - Uhh < h{h''^ + e''^ + c-^h + -t= + !^){\\u\\H^^r) + ll/llL^(r)). (3.57) 

Proof. The triangle inequality yields 

\\u^-uhh < \W - {Ihu'')\TA* + WhU^)\T^~Uhh. (3.58) 
The second term in the upper bound can be estimated using Lemmas 13. 2| 13.81 13.91 

\Whu'')\T,, - UhWl < ah{{Ihu'')\Ty^ - Uh, {Ihu'')\r^ - ut) 

= ah[[Ihu'')\T,, - u", [Ihu'')\T,^ - Uh) + Ohiu" - Uh, (/,iu'')|r^ - Uh) 

< h(h^n + ^1/2 + + + :^)||^||^.(r)||(/;,^-)|r, - uhh 

+ |a?i(u^ {Ihu'')\Y^ - Uh) - fh{{Ihu'')\r^ - Uh)\ 

< h{h'/^ + + cih + -jt= + :^)(||^,||^.(r) + ||/|U.(r))|l(/,X)|r, - «.|U. 
This results in 

||(/,u^)|r,-^i,.!l* < h{h^'^ + e''^ + c^h+^t= + l^){\\u\\H^^r) + \\f\\L^v)). (3.59) 

y e + c ye 

The error estimate (P37| follows from and 
□ 
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3.7. Further discussion. We comment on some aspects related to the main 
theorem. Concerning the analysis we note that the norm || • ||», which measures 
the error on the left-hand side of p.57p . is the standard SUPG norm as found in 
standard analyses of planar streamline-diffusion finite element methods. The analysis 
in this paper contains new ingredients compared to the planar case. To control the 
geometric errors (approximation of F by Th) we derived a consistency error bound in 
Lemma l3.9l To derive a continuity result (Lemma l3.8|) . as in the planar case, we apply 
partial integration to the term J^^{w^ ■ Vrft0)w/ids, cf. p.37p . However, different 
from the planar case, this results in jumps across the edges E € Sh which have to be 
controlled, cf. p.38p . For this the new results in the Lemmas 13.51 and 13.61 are derived. 

4 

These jump terms across the edges cause the term ^ in the error bound in p.57p . 

Consider the error reduction factor /i'^/^ + e^/^/i + c^h'^ + — ^ — h on the right- 
hand side of p.57|) . The first three terms of it are typical for the error analysis of 
planar SUPG finite element methods for PI elements. In the standard literature for 
the planar case, cf. [3], one typically only considers the case c > 0. Our analysis also 
applies to the case c = 0, cf. (|3.28|) . Furthermore, the estimates are uniform w.r.t. 
the size of the parameter c. For a fixed c > and s h the first four terms can be 
estimated by < /i'^/^, a bound similar to the standard one for the planar case. The 
only "suboptimal" term is the last one, which is caused by (our analysis of) the jump 
terms. Note, however, that < h^^'^ if <e holds, which is a very mild condition. 

The norm |j • ||* provides a robust control of streamline derivatives of the solution. 
Cross-wind oscillations, however, are not completely suppressed. It is well known that 
nonlinear stabilization methods can be used to get control over cross-wind derivatives 
as well. Extending such methods to surface PDEs is not within the scope of the 
present paper. 

The error estimates in this paper are in terms of the maximum mesh size over 
tetrahedra in w^, denoted by h. In practice, the stabilization parameter 6t is based on 
the local Peclet number and the stabilization is switched off or reduced in the regions, 
where the mesh is "sufficiently fine". To prove error estimates accounting for local 
smoothness of the solution u and the local mesh size (as available for planar SUPG 
FE method), one needs local interpolation properties of finite element spaces, instead 
of (|3.20p . (|3.2ip . Since our finite element space is based on traces of piecewise linear 
functions, such local estimates are not immediately available. The extension of our 
analysis to this non-quasi-uniform case is left for future research. 

Finally, we remark on the case of a varying reaction term coefficient c. If the 
coefficient c in the third term of (|2.4p varies, the above analysis is valid with minor 
modifications. We briefly explain these modifications. The stabilization parameter 5t 
should be based on elementwise values ct = maxxgTc(x). For the well-possedness of 
p.4p , it is sufficient to assume c to be strictly positive on a subset of F with positive 
measure: 

A := meas{x G F : c(x) > cq} > 0, 

with some cq > 0. If this is satisfied, the Friedrich's type inequality [12] (see, also 
Lemma 3.1 in [5]) 

\\v\\l.(r) < Ci.(|lVrz;|li.(r) + \\Vcv\\h^r)) ior all « G F 

holds, with a constant Cp depending on cq and A. Using this, all arguments in the 
analysis can be generalized to the case of a varying c(x) with obvious modifications. 
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With c„jiij := essinfxgrc(x) and Cmax •= esssupxgrc(x), the final error estimate takes 
the form 

ll""^ - uh\U < h{h'/^ + e'/' + cL.h + ^ + -^)(||u||^.(r) + ||/|U.(r)). 

4. Numerical experiments. In this section we show resuhs of a few numerical 
experiments which illustrate the performance of the method. 

Example 1. The stationary problem (|2.4p is solved on the unit sphere T, with 
the velocity field 



w(x) = (-2-2^1 - xi, a-i^l - xi, Of , 
which is tangential to the sphere. We set e — 10~^, c = 1 and consider the solution 

XiX2 ( ^3 

it(x) = arctan —= 

Note that u has a sharp internal layer along the equator of the sphere. The corre- 
sponding right-hand side function / is given by 



8e^/^(2 + £ + 2a;o)a;ia;2a;3 &£XiX2 + x/x^ + Xn(x^ — xV) f x^ 

7 1 — arctan 

7r(e 4-4x1)2 tt \y/e 



We consider the standard (unstabilized) finite element method in (|2.11l) and the 
stabilized method (|2.12p . A sequence of meshes was obtained by the gradual re- 
finement of the outer triangulation. The induced surface finite element spaces have 
dimensions N = 448, 1864, 7552, 30412. The resulting algebraic systems are solved 
by a direct sparse solver. Finite element errors are computed outside the layer: The 
variation of the quantities err]^2 — \\u — Uh\\L^(D)^ errui — \\u — and 
erri^inf = \\u — with Z) = {x £ F : jxaj > 0.3}, are shown in Figure l4Jl 

The results clearly show that the stabilized method performs much better than the 
standard one. The results for the stabilized method indicate a 0{h?) convergence 
in the L^{D)-noYm and L°°(£')-norm. In the iJ^(£')-norm a first order convergence 
is observed. Note that the analysis predicts (only) 0{h2) convergence order in the 
(global) L^-noim. 




2.5 3 3.5 4 4.5 2.5 3 3.5 4 4.5 2.5 3 3.5 4 4.5 

log,„(N) loaJN) log,„(N) 



Fig. 4.1. Discretization errors for Examvle\l\ 
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Figure 14.21 shows the computed solutions with the two methods. Since the layer 
is unresolved, the finite element method (|2.1H) produces globally oscillating solu- 
tion. The stabilized method gives a much better approximation, although the layer 
is slightly smeared, as is typical for the SUPG method. 





Fig. 4.2. Examvle U\ solutions using the stabilized method and the standard method 112.111 1. 

Example 2. Now we consider the stationary problem (|2.4p with c = 0. The 
problem is posed on the unit sphere F, with this same velocity field w as in Example[TJ 
We set £ = 10~^, and consider the solution 

X1X2 ( 

M(xj = arctan 



The corresponding right-hand side function / is now given by 



9,e^/^{2 + e + 2xl)xiX2Xz %exiX2 + \/xl: + xlix\ - xV) , / 2:3 
1 — arctan 



Note that for c = one looses explicit control of the L2 norm in || • Thus we 
consider the streamline diffusion error : 

errsD = ||wr • Vr(M - 

Results for this error quantity and for erri^inf ~ \\u — Uh\\L'^(D) a-re shown in Fig- 
ure 14.31 We observe a 0{h) behavior for the streamline diffusion error, which is 
consistent with our theoretical analysis. The L°°-norm of the error also shows a first 
order of convergence. 

Example 3. We show how this stabilization can be applied to a time dependent 
problem and illustrate its stabilizing effect. We consider a non-stationary problem 
()2.3p posed on the torus 



r = {(xi,X2,X3) I i^xl+xl-lf+xl = ^}. (4.1) 



We set e = 10 ^ and define the advection field 



V ^1 + ^2 
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streamline diffusion errors 




3.5 



Fig. 4.3. Discretization errors for Examvle [2| 



which is divergence free and satisfies w • nr = 0. The initial condition is 

Mo(x) 



X1X2 ^ I X3 

-arctan 



The function uq possesses an internal layer, as shown in Figure HT^ a). 

The stabilized spatial semi-discretization of (|2.3p reads: determine Uh = Uh{t) € 



such that 



with 



m{dtUh, Vh) + ah{uh, Vh) = for all Vh e Vf, 



(4.2) 



m{dtu,v) / dtuvds+ St / i9tu(w'^ • Vr^i') ds, 
ah{u,v) i Vr,,u • Vr,,wds + i / (w" • Vr^u)i; ds - / (w" • Vr^w)" ds 



_ ?i + w'^ • Vr. fi)w'^ • Vr, v ds. 



Note that a/i(-, •) is the same as a/i(-, •) in (j2.13p with c = 0. The resulting system of 
ordinary differential equations is discretized in time by the Crank-Nicolson scheme. 

For e = the exact solution is the transport of uo(x) by a rotation around the 
axis. Thus the inner layer remains the same for all t > Q. For e = 10~^, the exact 
solution is similar, unless t is large enough for dissipation to play a noticeable role. 
The space is constructed in the same way as in the previous examples. The spatial 
discretization has 5638 degrees of freedom. The fully discrete problem is obtained by 
combining the SUPG method in (|4.2p and the Crank-Nicolson method with time step 
5t = 0.1. The evolution of the solution is illustrated in Figure demonstrating a 
smoothly 'rotated' pattern. 

We repeated this experiment with i5t = in the bilinear forms m(-, •) and a/i(-, •) 
in (|4.2p . i.e. the method without stabilization. As expected, we obtain (on the same 
grid) much less smooth discrete solutions (Figure 14. 5p . 
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(c) (d) 
Fig. 4.4. Example{^ solutions for t = 0, 0.6, 1.2, 1.8 using the SUPG stabilized FEM. 




(a) (b) 
Fig. 4.5. Example\3i solutions for t = 1.2, 1.8 using the standard FEM. 



Remark 4.1. With respect to mass conservation of the scheme we note 
fohowing. For w/i = 1 in (|4.2|) we get, with Mhit) := /p^ uii{x,t) ds, 
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Here we used estimates from Lemmas 3.1 and 3.5. Using Lemma 3.6 we get 



E 



\uh\As<h--^( ^ 



EeSh 



u?ds 



< h ^{\\uh\\L^Tn) +H'^rnUh\\L^T,,))- 



Assume that for the discrete solution we have a bound ||M;i||i2(p^) + /ij| Vrftit/illLSj-p^) < 
c with c independent of h. Then we obtain {t)\ < h and thus \Mh{t) - M/i(0)| < 
ct/i, with a constant c independent of h and t. Hence, with respect to mass conserva- 
tion we have an error that is (only) first order in h. Concerning mass conservation it 
would be better to use a discretization in which in the discrete bilinear form in (|2.13p 
one replaces 



/ (w" • Vr;.it)wds - / (vif" • Vr;,u)wds by -/ (w^ • Vr^^w)" ds. (4.3) 



This method results in optimal mass conservation: -^Mhit) = 0. It turns out, how- 
ever, that (with our approach) it is more difficult to analyze. In particular, it is not 
clear how to derive a satisfactory coercivity bound. 

In numerical experiments we observed that the behavior of the two methods 
(i.e, with the two variants given in (|4.3p ) is very similar. In particular, the mass 
conservation error bound \M{t) — ^^{0)\ < cth for the first method seems to be 
too pessimistic (in many cases). To illustrate this, we show results for the problem 
described above, but with initial condition 



wo(x) 



1 H — arctan 



X3 



Mods = TT^ w 9.8696. 



Figure shows the quantity Mh{t) for several mesh sizes h. For t = we have, due 
to interpolation of the initial condition uq, a difference between Af/j(0) and /p Mods 
that is of order K^. For t > we see, except for the very coarse mesh with /i = 1/4 a 
very good mass conservation. 

Example 4. As a final illustration we show results for the non-stationary problem 
but now on a surface with a "less regular" shape. We take the surface given in 



{(xi,a;2,a;3) I (a;i - 2:3) +X2+x^^l). 



(4.4) 



We set £ ~ 10^^ and define the advection field as the F-tangential part of w = 
(—1,0,0)"^. This velocity field does not satisfy the divergence free condition. The 
initial condition is taken as uo(x) = 1. We apply the same method as in Example |3l 
The mesh size is 1/8 and the time step is 5t = 0.1. Figure H771 shows the solution for 
several t values. We observe that as time evolves mass is transported from the two 
poles on the right to the left pole, just as expected. Our discretization yields for this 
strongly convection-dominated transport problem a qualitatively good discrete result, 
even with a (very) low grid resolution. 
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Fig. 4.6. Total mass variation for Examvle l3\ 




(c) (d) 
Fig. 4.7. Example^ solutions for t = 0,0.5, 1.0,2.0 using the SUPG stabilized FEM. 
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